Density-functional embedding using a plane-wave basis 
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The constrained electron density method of embedding a Kohn-Sham system in a substrate system 
(first described by P. Cortona, Phys. Rev. B 44, 8454 (1991) and T.A. Wesolowski and A. Warshel, 
J. Phys. Chem 97, 8050 (1993)) is applied with a plane-wave basis and both local and non- 
local pseudopotentials. This method divides the electron density of the system into substrate and 
embedded electron densities, the sum of which is the electron density of the system of interest. 
Coupling between the substrate and embedded systems is achieved via approximate kinetic energy 
functionals. Bulk aluminium is examined as a test case for which there is a strong interaction 
between the substrate and embedded systems. A number of approximations to the kinetic-energy 
functional, both semi-local and non-local, are investigated. It is found that Kohn-Sham results can 
be well reproduced using a non-local kinetic energy functional, with the total energy accurate to 
better than 0.1 eV per atom and good agreement between the electron densities. 

PACS numbers: 71.15.Mb,71.15.Ap,71.15.-m 
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I. INTRODUCTION 

For the past two decades Density Functional Theory 
(DFT) [l|, @, E| has been one of the most powerful tools 
for the ab initio calculation of the physical and chemical 
properties of materials. Methods based on DFT make 
efficient use of computational resources, hence can gen- 
erally deal with larger and more complex systems than 
other ab initio methods. They also provide a simple in- 
terpretation of much of the many-electron physics of ma- 
terials in terms of ideas based on the electron gas. A 
number of implementations of DFT exist, which essen- 
tially differ in the approach taken to approximating the 
unknown density functional that describes the contribu- 
tion to the electronic energy that is not due to the exter- 
nal potential. 

The most successful of these methods is the approach 
first derived by Kohn and Sham [4], which uses DFT to 
identify the interacting electrons with a non-interacting 
electron gas, and then solves for the non-interacting elec- 
tron system. Most of the energy of the system is eval- 
uated exactly, with only a relatively small exchange- 
correlation contribution requiring approximation. In ad- 
dition, the exchange-correlation part of the functional 
can be well approximated by a simple analytic form, 
hence this has become the workhorse of accurate DFT 
Q. There is one main disadvantage that concerns us 
here. Within the standard Kohn-Sham approach the 
electron density is expressed as the electron density of 
a non-interacting, many-electron system, by obtaining 
the eigenstates of these non-interacting electrons; this 
requires 0(N 3 ) operations where N is the number of 
electrons present in the system. It is this scaling be- 
haviour that limits the size of system that the Kohn- 
Sham method can be applied to (currently less than 
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around 1000 atoms). 

A more direct and computationally cheaper approach 
is to minimise the total energy functional with respect 
to variations in the electron density, p(r) (eg Wang et al 
Q). In this form the cost of finding the minimum of the 
total energy does not depend on N provided the non- 
interacting kinetic energy functional T s [p] 3 is available 
as an explicit functional of the electron density. Unfortu- 
nately, this functional is not known, hence a direct min- 
imisation procedure must employ approximate forms of 
the kinetic energy functional as well as the exchange- 
correlation energy. Approximate expression are avail- 
able, but as T s [p] is generally an order of magnitude 
greater than the exchange-correlation energy they are 
generally not of sufficient accuracy for structural opti- 
misation, let alone chemical calculations. Another defi- 
ciency is that there is no obvious way of applying the 
familiar non-local pseudopotentials of Kohn-Sham meth- 
ods [5[ to these direct minimisation methods, although 
indirect methods have been proposed by Watson et al 
0, Anata and Madden [f| and Shah et al @|. 

Other methods have been investigated, such as a path 
integral formulation of Kohn-Sham theory [lfj| , and sev- 
eral approaches which formulate DFT as a 1 st order re- 
duced density matrix theory 0, [H| . The latter formula- 
tion takes advantage of the 'nearsightedness' (see Kohn 
[l2l ]) of the density matrix to solve for the ground state 
energy as an O(N) problem, with the non-interacting ki- 
netic energ y fu nctional evaluated exactly (eg Baroni and 
Giannozzi |i~3| . Hernandez et al [l4T|. Ordejon [l5j . and 
a recent review by Goedecker [161]). Although ab initio 
O(N) approaches of this form are successful, they are 
currently limited to insulators as the density matrix is 
long-ranged for systems with no band gap. 

A middle ground between the Kohn-Sham method and 
direct energy functional minimisation can be found via 
'embedding methods', and this is the approach consid- 
ered in this paper. In many cases the system we are 
interested in may be divided into two regions, / and //. 
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Region II is largely the same as a more simple system 
that may easily be solved for, whereas region I is where 
the interesting physics occurs. An example would be a 
defect in a crystal - region 77 would be a bulk crystal, and 
region I a small volume surrounding the defect. There is 
an obvious computational advantage in solving for region 
77 first, and then solving for region 7 taking into account 
the influence of region 77 in some way. This 'embedding' 
approach has received a great deal of attention, and a 
large number of methods have been presented in the lit- 
erature. No attempt is given to review all of these, but 
we refer to publications that describe general classes of 
methods. For the main part these differ in the space 
in which region 7 and region 77 are defined, fnglesfield 
[13] defined these regions in real space and constructed 
an exact embedding scheme that requires a knowledge of 
the Green function of the substrate (region 77), which is 
often prohibitively costly to calculate. Regions I and 77 
can also be defined more generally in Hilbert space, as 
described by Fisher UM- Another approach is discussed 
by Gutdeutsch et al where the 1 st order reduced 

density matrix is partitioned. All of these schemes are, 
to a greater or lesser extent, based on a wavefunction 
description of the embedding process. 

Here we apply a different procedure, first described by 
Cortona [2lJ and Wesolowski and Warshel 22], that car- 
ries out the embedding entirely at the density functional 
level. The essential idea is to express the total energy 
of the system in terms of two Kohn-Sham like systems 
with densities pi(r) and p2(r) (corresponding to regions 
/ and II), where the total electron density is given by 
p(r) = pi(r) + P2(r). This total energy is then min- 
imised by varying only the electron density pi(r), corre- 
sponding to fewer electrons than the whole system, and 
so a cheaper calculation. Potential applications of this 
method include defects in crystals, or adsorbates on sur- 
faces, since once a solution for the substrate is available 
the remainder of the calculation would only involve the 
Kohn-Sham representation of the electrons in the imme- 
diate vicinity of the defect or adsorbate. This partially 
frozen electron density method has been implemented by 
Cortona [Hj], Wesolowski et al [H, HI JM ill and (in a 
slightly different context) Govind et al 26] in a form that 
describes the coupling between the two subsystems via 
approximate kinetic energy functional. A localised set 
of basis functions was used in each case. 

In this paper this procedure is implemented within 
a plane- wave pseudopotential framework We in- 

vestigate the approach for metallic systems where the 
electrons in regions I and II are strongly interacting, 
whereas past applications have focused on insulating sys- 
tems and a relatively weak interaction. Our goal is to 
develop an embedding approach accurate enough, and 
computationally cheap enough, to aid the investigation 
of large scale defect and adsorbate systems. As a pre- 
liminary to this we test and analyse the approach by 
constructing a four atom cell of bulk fee aluminium by 
embedding one cubic sub-lattice of atoms within three 



others. Sections |TT] and IHII describe the implementation 
of the method, with approximate kinetic energy function- 
als described from the viewpoint of a plane- wave basis. 
Subsection IIV Al presents the results of the method ap- 
plied to bulk fee aluminium, with a number of different 
kinetic energy functionals. In subsection llV Bl a modified 
form of the method is applied in order to analyse the 
source of errors, again for different functionals. Another 
important issue in the direct application of DFT methods 
is the inclusion of non-local pseudopotentials. In subsec- 
tion IIV CI we present and justify the assumptions that 
must be made in order to employ non-local pseudopo- 
tentials within the embedding scheme. Rydberg atomic 
units are used throughout unless otherwise stated. 



II. PARTIALLY FROZEN ELECTRON DENSITY 

We begin with the familiar Hohenberg-Kohn total en- 
ergy functional E[p], expressed in the Kohn-Sham form 
1 

E[p\ = T s [p\ + J[p] + E xc [p] + J V ext (r)p(r)d 3 r (I) 

where T s [p] is the non-interacting kinetic energy func- 
tional, J[p] is the Hartree energy, V ext is the local ex- 
ternal potential (for the systems considered here it is 
given by the sum of local pseudopotentials at each atomic 
site) and E xc [p] is the exchange-correlation energy, in- 
cluding the kinetic-correlation contribution. In the stan- 
dard Kohn-Sham methodology the functional derivative 
of Eq. JT]) is taken with the total number of electrons 
constrained to be constant. By setting this equal to zero 
an Euler-Lagrange equation is obtained, which is iden- 
tified with the Euler-Lagrange equation for a system of 
non-interacting electrons in a specific external potential. 
Solving for this 'reference' system to yield the same elec- 
tron density as the interacting system is the essence of the 
Kohn-Sham implementation of density functional theory, 
and results in the electron density that gives the correct 
minimum of Eq. (fTJ , the ground state energy. 

The partially frozen electron density method breaks 
down this same functional into two non-interacting gases. 
The electron density of the system, p(r), is divided into 
two components so that p(r) = pi(r) + P2(r). One of 
these components (in what follows, /?2(r)) is taken to 
represent a part of the system that is expected to change 
very little (this will be qualified further on); the sub- 
strate. This is kept constant, and in what follows is ob- 
tained from a Kohn-Sham calculation so the kinetic en- 
ergy 7 s [p2] is known accurately. The total energy func- 
tional is then expressed in terms of the kinetic energy 
functional of pi(r) and an 'embedding kinetic energy' 
term, T™ add [pi, p 2 ] + T s [p 2 ] that takes into account the 
influence of the rest of the system. This gives the total 
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energy as [2l|, [22j 

E[ P ]= T s [ Pl ]+Tr dd [pi,P2}+T s [p 2 ] 

+ J[p] + E xc [p\ + [ V ext (r)p(r)d 3 r (2) 



where the non- additive part of T B [p], T™ add [pi, p 2 ], is de- 
fined as 

TT dd [pi,P2] = T s [pi + p 2 ] - T.fa] - T s [p 2 }. (3) 

Minimising Eq. ([2]) with respect to variations in P \ (r) 
only, with the substrate density p 2 (r) constant, and a 
constraint of constant total number of electrons in pi (r) , 
results in the Euler-Lagrange equation 
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+ V KS [p;r] =h (4) 



where p is an arbitrary constant reflecting the fact that 
for a fixed number of electrons the functional derivative is 
defined to within an additive constant only. In the same 
manner as for the Kohn-Sham case, this leads to the 
pi(r) being the solution of the 'Kohn-Sham' equations 
associated with Eq. ((4]) at self consistency, but with an 
effective potential given by 

dpi 

= %[p;r] + f^-«\5) 



8 Pi 



Spi 



If an exact expression for the kinetic energy functional 
was available this prescription would provide a ground 
state energy and electron density exactly equivalent to 
the Kohn-Sham scheme for the total system, with one 
additional limitation. Since pi (r) takes the form 



Pl (r)=^^|^(r)| 5 



(0) 



it is positive, hence the trial densities that are searched 
to minimise the total energy in Eq. (fT|) satisfy p(r) > 
p 2 (r), and the true ground state energy and density is 
obtained only if the ground state density satisfies this 
inequality. In practise p 2 (r) is chosen such that this is 
not a significant restriction, and this is true for the test 
cases considered in section ITVl This constraint may also 
be relaxed by applying a 'Freeze and Thaw' procedure, 
as discussed at the end of subsection IIVBI 

Since p 2 (r) is taken as already known, this method of 
obtaining the electronic structure need only solve for the 
electrons present in pi(r), a smaller number (in many 
cases considerably smaller) than is present in the entire 
system. However, in order to apply this approach the 
term T™ add {pi, p 2 ], or non-additive kinetic energy f27\,\2§^, 
and its functional derivative in Eq. @ and Eq. ([5]) are 
required. Since no explicit form is available approximate 
kinetic energy functionals are employed in Eq. ([3]) to 



provide an approximate non-additive kinetic energy func- 
tional. To clarify, we approximate all of the functionals 
on the RHS of Eq. whereas in Eq. @ the functionals 
T s [pi] and T„[p2] are exact. This approximation to the 
non-additive kinetic energy is the only additional source 
of error introduced by the method, but T™ add is expected 
to be far smaller than the total kinetic energy (it is zero 
if pi(r) and p 2 (r) do not overlap) for most reasonable 
divisions of the electron density, and it could be hoped 
that some error cancellation will occur. 

The non-additive kinetic energy has been investigated 
as a test of the quality of a number of kinetic energy func- 
tionals by Lacks and Gordon (27[, who compared T™ add 
for Helium and Neon calculated from approximate func- 
tionals with Hartree-Fock results. They conclude that 
the fractional error of T™ add is greater than for the kinetic 
energy of the whole system, implying that for weakly 
interacting subsystems the error cancellation is limited. 
Whether this is the case for more strongly interacting 
subsystems and for the functionals applied here will have 
a direct influence on the accuracy of our results. In ad- 
dition it should be remembered that even if the approx- 
imate functional gives the correct kinetic energy for the 
true density, this may not be a minimum of the total 
energy with respect to variations in the density. 



III. APPROXIMATE KINETIC ENERGY 
FUNCTIONALS 

Approximations to the kinetic energy functionals are 
used to construct T nadd \p 1 , p 2 ] in Eqs. &J5]). Many are 
available in the literature (eg Thakkar 12911 . Wang et al 
[H, Garcfa-Gonzalez et al [33], Herring |3l|), and these 
have been assessed in a number of environments ranging 
from isolated atoms to bulk systems to molecular inter- 
actions and surfaces. The majority of these assessments 
investigate the ability of the functionals to reproduce ac- 
curate kinetic energies from accurate electron densities 
calculated by other means. A smaller number of stud- 
ies have examined the ability of approximate functionals 
to produce accurate electron densities and energies when 
the total energy is minimised using the functionals them- 
selves, and little work has been published on the success 
of these approximations in reproducing accurate func- 
tional derivatives. We have therefore chosen to examine 
a range of functionals. 

Although an analytic gradient expansion exists for 
T s [p] , convergence cannot be achieved for systems where 
the density decays exponentially 2]. This has been at- 
tributed to the expansion taking the form of an asymp- 
totic series, as described by Pearson and Gordon [32j |. 
although this has not been shown analytically. In order 
to overcome this difficulty in improving the local density 
approximation (LDA) to T s [p] , many authors have taken 
a similar approach to the Generalised Gradient Approxi- 
mation (GGA) to the exchange-correlation functional, by 
carrying out a partial re-summation via an enhancement 
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factor 29]. The kinetic energy is approximated by 



T?PP[p] = |(37T 2 )t / p*F(t)d 3 r 



where 



t = 



(7) 



(8) 



hence the kinetic energy is expressed as a semi-local func- 
tional of the electron density and its gradient. The func- 
tional derivative of a general semi- local density functional 



G[ P ] 



)d 3 



is given 



by 1 

SG 



dg 
dp 



dg 
'dVp 



< dg 
dV 2 P 



(9) 



(10) 



which may be truncated at second order for the func- 
tional given in Eq. ([7]), since the kernel is a function of 
p and Vp only. 

The direct application of Eq. (fit))) to the semi-local 
functional, Eq. (7J, yields 
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an unwieldy expression that involves highly non-linear 
terms. With a plane-wave basis difficulties arise due to 
p(r) being defined on a real space grid, with gradients 
conventionally obtained via the Fast Fourier Transform 
(FFT). Since Eq. (fTTj) is non-polynomial in p(r) it cannot 
be represented by a finite Fourier space, and using any 
finite space results in aliasing errors. For example, the 
term Vp.V|V,o| results in particularly large errors, and if 
Eq. (fTTj) is applied directly this causes errors in the re- 
sulting potential to propagate through further iterations, 
preventing convergence. 

Exactly this problem manifests itself in the applica- 
tion of the GGA with a plane- wave basis, as discussed 
by White and Bird[33j. In this case convergence is also 
affected, though not as severely due to the exchange- 
correlation energy being an order of magnitude smaller 
than the kinetic energy. We follow the same approach 
as White and Bird to find a more stable expression for 
the kinetic energy functional derivative. The functional 
is discretised as a sum of contributions at each real space 
grid point, 



T app [p] = -(3tt 2 V 



N ^ 

R 



(12) 



where f2 is the volume of the unit cell, N is the number of 
grid points, and the subscript R denotes the quantity at 
the grid point R. The 'functional derivative' (in fact a to- 
tal derivative) can then be written in terms of the partial 
derivative of Eq. (| 1 2[) with respect to /?r and VpR, which 
are considered as independent variables. This yields the 
expression 
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V. 
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which is analytically equivalent to Eq. (|11[) in the limit 
N — > oo. Expression (| 13[) for the functional derivative 
is applied in our calculations, with all gradients calcu- 
lated via the FFT, and is found to remove the instability 
inherent in Eq. (|lip . 

It is important to note that although Eq. (fT"3f defines 
the functional derivative in a manner consistent with a 
finite real space grid, it does not ensure that an accurate 
representation of the functional derivative is obtained for 
a given electron density. This is immediately apparent 
if an exponentially decaying electron density is consid- 
ered. In this case the second and third terms in Eq. (fT3"j) 

should result in the expressions I Vp £ l I and being 

Pr, PR- 
constant due to a cancellation of the exponentials. How- 
ever, since the numerator in both these terms will in fact 
be the gradient of a trigonometric representation of the 
electron density an exact cancellation will not occur, and 
for small values of pr errors in these terms (and the re- 
sulting functional derivative) will not be small. This is 
found to prevent convergence if the derivative of F(t) in 
Eq. (|T3")) is large, but has no effect for the functionals 
considered here. 

A number of enhancement factors available in the liter- 
ature [H, [29| were investigated (specifically the Thomas- 
Fermi approximation and 1 st order gradient expansion 
, and the functionals constructed by Thakkar [2i| , Vi- 
tos et al [HI , DePristo and Kress [35[ , Ou-Yang and 
Levy H, Lee et al H3, Perdew and Wang [H, |Ha S3] 
and Lembarki and Chermette [Ilj]), but their overall be- 
haviour was found to be similar. We therefore present 
results for two enhancement factors only. The first is 
that of Perdew and Wang, Fpwse [lil , 
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where 
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(14) 
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We choose this enhancement factor since Lacks and Gor- 
don (27j found it gave the best approximation to the non- 
additive kinetic energy of a number of related functionals. 
The second enhancement factor is 



Ftf- 



XvW 



(s) 



1 



(16) 
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which for A = is the Thomas-Fermi (TF) approxima- 
tion and for A — i the gradient expansion truncated at 
first order. The parameter A is taken as a free param- 
eter in order to optimise the results of the calculations 
(this is not a completely empirical approach and some 
theoretical justification is available 0, EH). 

A number of more general non-local approximations 
to the kinetic ener gy f unctional are available, such as lo- 
cal density scaling [43| and weighted density approxima- 
tions (eg Chacon et al [iij .Garcia-Gonzalez et al [30l]). 
Although these approaches are accurate they are com- 
putationally expensive, hence we choose a simpler form 
which has been found to accurately reproduce energies 
and densities for some bulk systems jfjj. This functional 
is from the family of approximations introduced by Wang 
and Teter [45J, Perrot [46j| and by Smargiassi and Mad- 
den [I?]], summarised and generalised by Wang et al @. 
These take the form 



rpl, 



~\p\ = -(3tt 2 )3 / pid 3 r- J p2V 2 p*d 3 r 

+ J p a w a {r-r')p a d 3 rd 3 r' (17) 



where the first term is the TF functional, the second term 
is the von Weizsacker (vW) functional and the third term 
is defined such that the entire functional has the correct 
linear response for a homogenous, non-interacting elec- 
tron gas. A number of different values have been pro- 
posed for the parameter a, each with its own justifica- 
tion as discussed by Wang et al @. For a plane- wave 
representation of p(r) Eq. (fTT|) takes the form 



-^nloc 



F 2 2 
19 P- 



(18) 



where the powers of the electron density are taken before 
transformation to reciprocal space, and is the volume 
of the unit cell. The linear response correction, w a (g) is 
given by 



k 2 



1 



3a 2 Po a_1 
1 



4a 2 p 2 Q a 
1 
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2a 2 p 2 {a ~ 1} XMnd(g) 



(19) 



with XLind the Lindhard susceptibility function in recip- 
rocal space, 
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Bird's treatment of the semi-local functionals. Introduc- 
ing an infinitesimal change in p(r) at a grid point results 
in 



8T?° c \p] 
5p 
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but 



where kf = (37rpo) 3 is the Fermi vector for the average 7T 
electron density po, and 77 = The functional deriva- 
tive can be obtained by writing Eq. (|17[) as a discrete 
sum (double sum for the linear response correction) over 
the real space grid, in the same manner as White and 



Equation (f21j) was found to be stable for a = 2 . 
became unstable for other values. This behaviour can be 
ascribed to the fact that in the limit of g — > Eq. (fTT)) 
gives the exact second order gradient expansion only for 
a = \ [6]. For all other values this not the case, hence 
the convergence problems associated with exponentially 
decaying, low electron density regions described earlier 
become significant. 



IV. EMBEDDING TESTED ON BULK FCC 
ALUMINIUM 

A. All terms of T™ add approximate 

To investigate this partially frozen density approach 
we examine fee aluminium with a 4 atom cubic unit cell. 
The 3 face-centred atoms are taken to be the substrate 
system, and this structure is solved to provide pi (r) and 
T s [p2]- A plane- wave basis pseudopotential approach is 
used, with a lattice constant of do = 4.05A, a plane- wave 
cut-off of 200 cV, 35k points in the irreducible wedge 
of the Brillouin zone, and the Goodwin-Needs-Heine [48| 
local pseudopotential (a non-local pseudopotential is con- 
sidered in subsection IIV C[) . Exchange-correlation is de- 
scribed by the LDA. 

Once this substrate is constructed the embedded 
Kohn-Sham calculation is carried out as a standard 
plane-wave basis calculation, but with the trial poten- 
tial given by Eq. ([5]) and the total energy given by Eq. 
@ • The non-additive kinetic energy in Eq. is given by 
Eq. with all of the terms on the RHS approximate. 
The unit cell, Brillouin zone sampling and other param- 
eters of the calculation are chosen to be the same as the 
substrate calculation. It should be made clear that for 
the substrate calculation we are solving for the lattice of 
3 face centred atoms and their accompanying electrons, 
but for the embedded calculation we are solving for the 
entire fee system, but only the electrons associated with 
the embedded (corner) atom are provided with a Kohn- 
Sham representation. 

We discuss results for 3 approximate functionals. 
Equation (fl"6f was applied with A = | (denoted 



tf-^vw)i smce this value was found to provide a useful 
compromise between accuracy of the value of the func- 
tional itself and the functional derivative. We also em- 
ploy Eq. (fl"4"|) {T PW86 ) and the non-local linear response 
corrected functional, Eq. (1 1 T[) . with a = h (denoted 
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Peak error/ 



Functional 


E/eV AE/eV R/% 


xl0~ 3 A -3 


T TF _4 vW 


-58.722 -0.392 4.206 


28.750 


Tpw86 


-58.052 0.277 6.390 


31.491 


rpnloc 
2 


-58.337 -0.008 4.614 


20.051 


Kohn Sham 


-58.329 





TABLE I: Total energy per atom E, and errors in energy 
AE and electron density. Results obtained with embedding 
scheme described in subsection IIV Al 

T? loc ). In Table U the total energy per atom of the em- 

2 

bedded calculations are given, together with the result of 
a full Kohn-Sham calculation carried out with the same 
basis and pseudopotcntial. T™'°° gives by far the best re- 

2 

suit for the total energy, although the extreme accuracy 
of this value is probably spurious since similar systems 
consistently provide errors of order ~ 0.1 eV/atom. The 
total energy for the other functionals is not as accurate, 
and results for other semi-local functionals that are not 
given here are similar or worse. 

Fig. QJi, Fig. [Tp and Fig. [TJ: show the electron density 
obtained using the T TF _± vW , Tpwm and T"' oc function- 

9 2 

als respectively. These figures show the electron density 
along a line in the [111] direction between two embed- 
ded atoms on the corners of the cubic unit cell. Table 
U gives the error in the electron density, compared with 
the Kohn-Sham results, over the whole unit cell. This 
is quantified as the peak error and the average absolute 
error, R, given by 

R = J |p(r) - p KS {v)\d z v/ J p KS (r)d 3 r, (22) 

where p KS (r) is the Kohn-Sham electron density. From 
these results it is apparent that the electron densities 
show the correct behaviour, reproducing the structure of 
the Kohn-Sham electron density reasonably well. How- 
ever, errors appear near the atomic sites of the embedded 
atoms as well as in in the region where pi(r) and /02(f) 
have the greatest overlap. Considering the errors in both 

the electron density and total energy, Tf loc provides the 

2 

most accurate reproduction of the Kohn-Sham results. 

B. One term of T™ dd approximate 

In this section we approximate the non-additive kinetic 
energy in a form that represents the kinetic energy of 
the entire system as an approximate kinetic energy func- 
tional. However, the densities px(r) and P2( r ) are still 
represented as Kohn-Sham systems. This corresponds 
to the many applications of approximate kinetic energy 
functionals to the direct minimisation of the total energy 
function (see Goedecker with no Kohn-Sham repre- 
sentation. 




1 2 3 4 5 6 7 

x(k) 




1 2 3 4 5 6 7 




1 2 3 4 5 6 7 



x(A) 



FIG. 1: Electron density in [111] direction for fee aluminium 
using the embedding scheme described in subsection IIV Al 
Embedded atoms are at 0.00 and 7.01 A. The dashed line 
shows the Kohn-Sham result and the solid line the embedding 
results for a) T TF _± vW functional, b) Tpwm functional and 

c) Tf oc functional (see text). 
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At first representing the electron density as the sum of 
two Kohn-Sham representations may seem like a waste 
of computational effort since the evaluation of the ap- 
proximate kinetic energy functional for the entire system 
does not require anything more than the electron density 
itself. However, it could be useful in two different ways. 

First it allows non-local pseudopotentials to be applied 
directly to each part of the system within the Kohn-Sham 
framework (see subsection II V Cp . Second the results will 
tell us whether the error in T™ add [pi, p 2 ] is greater or 
smaller than the error in T® pp [pi + p 2 \ ■ This second point 
is important since for the method described in the pre- 
vious section to be useful the error in the non-additive 
kinetic energy must be smaller than the error in the to- 
tal kinetic energy, as described by the approximate func- 
tionals. Some cancellation of errors must take place in 
T" add [pi, P2], and its functional derivative, for this to be 
the case. From the previous subsection it is apparent 
that cancellation occurs only to a limited degree, and a 
similar conclusion has been reached when addressing the 
accuracy of kinetic energy functionals when used to eval- 
uate interaction energies, both for gradient expansions 

28 1 and semi-local enhancement factor approximations 

27|. 

In order to obtain the required functional the non- 
additive kinetic energy is defined as 

Tr dd [pi, P2] = T^[ Pl + p 2 ] - T s [p x ] - T s [p 2 ], (23) 

where the first term on the RHS is an approximate func- 
tional, and the remaining terms are exact. In the pre- 
vious subsection all of the functionals in this expression 
were evaluated using approximate functionals. Equation 
@ is applied as before, with T s [/5i] and T s [p2\ exact, but 
with T nadd given by Eq. . The functional derivative 
of the non-additive kinetic energy becomes 



STI ladd 



px,p 2 ] 6T?*\pi + p2] ST s [ Pl ] 



Spi 



Spi 



Spi 



(24) 



where the first term on the RHS is approximate, and the 
second exact. To obtain the second term in Eq. (j2"4"|) 
the method of Bartolotti and Acharya [!§] is applied. 
They derive an expression for the functional derivative by 
replacing the self-consistent potential within the Kohn- 
Sham equations with the self consistent potential in 
terms of the associated Euler-Lagrange equation, 



■V 2 V„(k)- 



ST s [ Pl 
Spi 



^„(k)=e n (k)^„(k) (25) 



which gives 



ST S [pi 
Spi 



V 2 ^n(k) 
V-n(k) 



e„(k)+// 



(26) 



where p! is the associated Fermi energy. It should be 
noted that the Euler-Lagrange equation can be used di- 
rectly to obtain the functional derivative in terms of the 







Peak 


error/ 


Functional 


E/eV AE/eV R/% 


xlO" 


-3 A -3 


T TF _i vW 


-58.769 -0.440 3.922 




25.050 


TpwS6 


-59.680 -1.351 7.260 




57.064 


rpnloc 
2 


-58.411 -0.082 1.204 




9.180 


Kohn Sham 


-58.329 







TABLE II: Total energy per atom E, and errors in energy 
AE and electron density. Results obtained with embedding 
scheme described in subsection II V Bl 



trial potential, but this was found to cause convergence 
difficulties. To determine p! we use 

TM = J Pl (r)^ild 3 r, (27) 

as derived by Liu and Parr [50j . 

In effect this second approach corresponds to minimis- 
ing the energy functional of the entire system with the 
kinetic energy represented by the approximate kinetic 
energy functional, but with the additional constraint of 
p(r) > /°2(r)- Calculations were carried out within this 
scheme, with all parameters equivalent to those in the 
previous subsection. Table |TT] shows the total energy per 
atom, and the associated errors. Generally, the errors in 
the energy for the semi-local functionals are considerably 
worse than in Table []] as demonstrated by the results for 
Tpw86i but the error for T TF _± vW is similar to the first 
method. The non-local functional gives by far the most 
accurate results, with the energy accurate to better than 
0.1 eV. Fig. [2ji, Fig. Eb and Fig. Efc show the electron 
density resulting from these calculations, for T TF _± vW , 

Tpws6 and T± loc respectively. It is immediately apparent 

2 

that the non-local functional provides a far more accurate 

electron density than either of the semi-local functionals, 

and a more accurate electron density than the approach 

described in subsection IIVAI 

In subsection IIV Al errors in the total energy functional 

are introduced by the non-additive part of the kinetic 

energy, whereas in this section the errors are due to the 

approximate total kinetic energy. Bearing this in mind 

our results indicate that, for the non-local functional, 

T[pi + P2] is particularly well described by T" ioc , but 

2 

T[pi] and T[p 2 ] are not quite as accurate (these terms 
only appear in the non-additive kinetic energy). This 
difference in the accuracy of the same functional applied 
to different electron densities can be ascribed to the dif- 
ferent characters of the electron densities. Fig. [3] shows 
embedded, substrate and total electron densities result- 
ing from a calculation of the type described in this sec- 
tion, carried out with the non-local functional. Electron 
densities are shown along a line in the [110] direction be- 
tween two embedded atoms on opposite corners of one 
face of the cubic unit cell (this line is chosen as it in- 
cludes a substrate atom along its path). It is apparent 
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FIG. 2: Electron density in [111] direction for fee aluminium 
using the embedding scheme described in subsection II V Bl 
Embedded atoms are at 0.00 and 7.01 A. The dashed line 
shows the Kohn-Sham result and the solid line the embedding 
results for a) T TF _± vW functional, b) Tpwm functional and 

c) Tf oc functional (see text). 



that p\ is far from homogeneous, and falls close to zero 
near the sites of substrate atoms. This is expected to 
be accompanied by worse performance of T™ loc [pi] than 

2 

Tf [pi + P2], since this approximate functional is dc- 

2 

rived using arguments based on the linear response of a 
homogenous electron gas. It has previously been found 
that this functional is particularly successful for the total 
electron density of bulk aluminium with the pseudopo- 
tential applied here @. Similar considerations lead us 
to conclude that in general for semi-local functionals the 
difference T^p[ Pi + p 2 ] - T^p[ Pi ] - T^P[p 2 ] is more ac- 
curate than T s app [pi + p 2 ] (or about the same accuracy 
for the particular form T TF _± vW ). 

It is worth noting that the results obtained here for 
Aluminium show closer agreement with Kohn-Sham re- 
sults than those obtained by Wang et al Q using a 
slightly generalised form of the non-local functional, the 
same pseudopotential and direct minimisation of the to- 
tal energy functional. This is probably due to the ad- 
ditional constraint of p(r) > pa(r) present in the calcu- 
lation performed here preventing an over-relaxation of 
the electron density. This constraint can be relaxed by 
swapping the substrate and embedded systems after self- 
consistency has been reached, reaching self consistency 
again and repeating this until convergence of the total 
system is reached. Performing this 'Freeze and Thaw' 
[22I [H| procedure offers no particular advantage for the 
test system investigated in this paper and when imple- 
mented did not influence any of the conclusions given. 
The only significant consequence of performing a 'Freeze 
and Thaw' calculation was a 'drift' of the electron den- 
sities away from their associated atoms. For the non- 
local and T TF _4 vW functionals the two electron densi- 
ties remain localised on their host atoms, whereas for 
the Tpw86 (and other enhancement factor functionals) 
the two electron densities evolved to fill the entire unit 
cell with no association with the host atoms of region / 
or II apparent. 



C. Non-local pseudopotentials 

In past applications of the Hohenberg-Kohn theorem 
directly to the minimisation of the total energy functional 
the main problem addressed has been the inaccuracy of 
the available approximate kinetic energy functionals. An 
additional problem is how to make use of non-local pseu- 
dopotentials when no Kohn-Sham representation of the 
electron density is available. This problem arises due to 
the Hohenberg-Kohn theorem being strictly applicable 
only for a local external potential [HI], [HJ . 

For the embedding method applied here no Kohn- 
Sham representation of the total electron density is avail- 
able, so it is not immediately apparent how a non-local 
pseudopotential can be applied. This issue has been ad- 
dressed previously, for example Shah et al Q provide 
an ad hoc prescription, by taking the square root of the 
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FIG. 3: Electron density in [110] direction for fee aluminium 
using the embedding scheme described in subsection II V Bl and 
the T£ ioc functional (see text). Embedded atoms are at 0.00 

2 

and 5.73 A, and substrate atom is at 2.86 A. The dashed 
line shows the embedded electron density, pi, the dotted line 
the substrate electron density, p2, and the solid line the total 
electron density, p = pi + P2- 

density as the basic variable, while Watson et al [7| and 
Anata and Madden Q describe a procedure to obtain a 
new local pseudopotential from a non-local one. Here we 
make two assumptions. 

First, that the Kohn-Sham density matrix 7 s (r,r'), 
given by 

7s (r,r')=^^V l -(rM(r') ! (28) 

i 

(where u>i is the occupation number for state i) is an 
accurate enough representation of the actual many-body 
1 st order reduced density matrix (this is the assumption 
implicit in the normal application of non-local pseudopo- 
tentials within the standard Kohn-Sham scheme). As- 
suming this to be a valid approximation the external po- 
tential energy is given by 

E ext [ ls {r,r')] = J t>(r',r) 7s (r,r')d 3 r'd 3 r (29) 

where w(r',r) is an external non-local potential. Second 
we take advantage of the fact that the electron density 
component pi (r) is mostly localised near the atomic sites 
in the substrate, and p\(y) is expected to remain localised 
around the site of the embedded atoms. Bearing this in 
mind we take the exact expression 

p(r)= /9l (r)+p 2 (r) (30) 

and generalise it to the approximation 

7 s (r,r') w 7sil (r,r') +7 s , 2 (r,r'), (31) 



Peak error/ 



Functional 


E/oV AE/eV R/% 


xl0~ 3 A -3 


T TF _i vW 


-57.148 -0.120 4.349 


25.267 


TpwS6 


-56.262 0.767 7.137 


41.953 


rpnloc 

2 


-56.940 0.088 6.343 


43.735 


Kohn Sham 


-57.028 





TABLE III: Total energy per atom E, and errors in energy 
AE and electron density. Results obtained with embedding 
scheme described in subsection IIV Al using a non-local pseu- 
dopotential. 



where 7^1 and 7si 2 are the 1 st order reduced Kohn-Sham 
density matrices corresponding to electron densities p\ (r) 
and P2(r). Of course this can only be exact, over all 
space, if there is no overlap between the electron densi- 
ties, as can easily be deduced from the fact that there is a 
non-additive component to the kinetic energy functional 
in the first place. 

Approximation (|31| immediately leads to 



E ex t[ls\ ~ J Vi oc (r)p(r)d 3 r + Vniocfoss] + V n i oc [73,2] 

(32) 

where the first term is the local part of the potential, 
and the second and third terms are the contributions of 
the embedded and substrate systems due to the non-local 
part of the pseudopotentials of all the atoms. Conven- 
tional norm-conserving pseudopotentials are non-local 
only between points on the surface of spheres centred on 
each atomic site, and only for spheres with a radius less 
than a certain value, r c . This implies that the approxi- 
mation in Eq. (f3"Tj) need only be accurate between points 
on each such sphere. Since /92(f) is largely localised near 
the atomic sites in the substrate, and pi(r) is expected 
to remain localised near the embedded atomic sites, Eq. 
(|32[) provides a reasonable approximation. 

Calculations were carried out for fee Al with all param- 
eters as before, and using the approach given in subsec- 
tion IIV Al but with a Kerker non-local pseudopotential 
[53J. Fig. 2] and Table IHII show results for a calcula- 
tion with p2(r) frozen and for the T TF _± vW , Tpwa6 and 
non-local functional. These results should be compared 
with those of Fig. [1] and Table Q] For the total ener- 
gies both T TF _4 vW and T" Zoc yield results essentially as 
9 2 

accurate as those obtained for the local pseudopotential 
in subsection IIV Al while Tpwse is less accurate. Fig. 0] 
shows a higher peak error in the electron densities, but 
it is difficult to attribute this to a particular aspect of 
the approximations inherent in the implementation of the 
non-local pseudopotential. We conclude that a non-local 
pseudopotential can be used within this method with no 
significant loss of accuracy. 
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FIG. 4: Electron density in [111] direction for fee aluminium 
using the embedding scheme described in subsection llV Al and 
a non-local pseudopotential. Embedded atoms are at 0.00 and 
7.01 A. The dashed line shows the Kohn-Sham result and the 
solid line the embedding results for a) T TF _± vW functional, 

b) Tpwm functional and c) T? loc functional (see text). 



V. DISCUSSION AND CONCLUSION 

Wc have implemented the partially frozen electron 
density approach of Cortona [21| and Wesolowski and 
Warshel [23 . |23| using a plane-wave basis, both local 
and non-local pseudopotentials, and for a metallic sys- 
tem. Although some numerical instabilities are intro- 
duced by using a plane- wave representation these are con- 
trolled using methods previously developed for exchange- 
correlation energies and potentials. Several approxi- 
mations for the kinetic energy functional are consid- 
ered, including the semi-local enhancement factor ap- 
proximations that have previously been applied within 
this method. In addition a modified Thomas-Fermi/von 
Weizsacker functional and a non-local functional are im- 
plemented. 

A lattice of aluminium atoms is embedded into a sub- 
strate lattice of atoms to create bulk fee aluminium, and 
we find that the semi-local functionals result in total en- 
ergies within ~ 0.2 — 0.5 eV per atom of the Kohn-Sham 
result, and the non-local functional results in an energy 
differing by < 0.1 eV per atom. Kohn-Sham electron 
densities are reproduced reasonably closely. The non- 
local functional performs best with a peak error of ~ 20 
milli-electrons A -3 while the semi-local functionals result 
in errors of ~ 30 milli-electrons A -3 . Calculations per- 
formed with a non-local pseudopotential produce results 
for the total energy and electron density which are of 
comparable accuracy to the local pseudopotential case. 

Calculations are also carried out for an altered form of 
the method where part of the kinetic contribution to the 
embedding potential and energy is obtained exactly. This 
corresponded to performing a Hohenberg-Kohn minimi- 
sation of the total energy expressed in terms of the elec- 
tron density with an approximate kinetic energy func- 
tional applied to the entire system. As implemented here 
this minimisation allows the use of non-local pseudopo- 
tentials, and introduces the constraint p > P2 where /?2(r) 
is some reference substrate system. Results are worse 
than for the true embedding scheme for all the semi-local 
functionals, with the exception of T TF _± vW which gives 
a slightly greater error in the energy, but a slightly im- 
proved electron density. From this it seems reasonable to 
conclude that the approximate non-additive kinetic en- 
ergy is more accurate than the approximate total kinetic 
energy for these functionals. The non-local functional 
gives the most accurate results, with the energy accurate 
to < 0.1 eV atom -1 and the electron density accurate 
to < 10 milli-electrons A -3 . This suggests that the non- 
local functional produces the most accurate representa- 
tion of the both the value of the kinetic energy functional 
and its functional derivative, and that the errors in the 
functional derivative are greatest when the electron den- 
sity is low. 

We interpret the extremely good agreement between 
the Kohn-Sham electron density and embedding results 
found in subsection llV Bl ( see Fig.^fc) in comparison with 
that found in subsection IIV Al (Fig. [Th) for the non-local 
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functional to be due to the success of this functional in 
describing bulk aluminium as discussed by Wang et al 
Q. For systems where the total electron density is far 
from homogeneous (eg surface/adsorbate) this success of 
the method of subsection IIVBI is not expected to hold. 
In future applications the non-local functional and the 
method of subsection IIV Al are expected to provide the 
most accurate reconstruction of the full Kohn-Sham re- 
sult. 
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